Overview

This markdown includes data, code, figures, and statstics for the physiology analysis of the Porites Rim Bleaching experiment.

Air brushing Protocol

Symbiont density Protocol

Chlorophyll Protocol

Total Protein Protocol

Carbohydrate Protocol

Loading libraries and metadata files

# load packages
library(dplyr)
library(tidyverse)
library(readr)
library(stringr)
library(gridExtra)
library(grid)
library(ggplot2)
library(lattice)
library(Rmisc)
library(ggpubr)
library(lsmeans)
library("reshape")
library("arsenal")
library(RColorBrewer)
library(lme4)
library(lmerTest)
library(car)
library(effects)
library(ggfortify)
library(cowplot)
library(vegan)
library(corrr)
library(ggcorrplot)
library(GGally)
library(broom)
library(cowplot)
library(RVAideMemoire)
library(mixOmics)
library(factoextra)
library("segmented")
library("plotrix")
library("lubridate")
library(lsmeans)

# Load metadata
vial.meta <- read.csv("../data/Physiology/Vial_Metadata.csv")
frag.meta <- read.csv("../data/Physiology/Airbrush_Metadata.csv")

vial.meta$coral.tp <- paste(vial.meta$Fragment.ID, vial.meta$Timepoint, sep = "-")
frag.meta$coral.tp <- paste(frag.meta$Fragment.ID, frag.meta$Timepoint, sep = "-")

master.meta <- merge(vial.meta, frag.meta, by = "coral.tp")

Endosymbiont Density

#Import Data
Zoox <- read.csv("../data/Physiology/Symbiont_Counts/Symbiont_Counts.csv")

#Calculating cells/Larvae from cells/count
Zoox$Average <- rowMeans(Zoox[ ,3:8]) #averaging counts
Zoox$Cells.mL <- (Zoox$Average/Zoox$Num.Squares)/0.0001 #number of cells per mL - volume of haemocytometer

# Merging metadata

zoox.meta <- merge(Zoox, master.meta, by = "Vial")

#Accounting for homogenate volume 
zoox.meta$Cells.mL.vol <- zoox.meta$Cells.mL * zoox.meta$Homogenate_Vol.mL

#Normalizing to fragment surface area
zoox.meta$Cells.cm2 <- zoox.meta$Cells.mL.vol/zoox.meta$Surface_Area_cm2
zoox.meta$Cells.cm2.x6 <- zoox.meta$Cells.cm2 / 1000000

#Removing metadata columns
zoox.final <- zoox.meta %>%
  dplyr::select("Fragment.ID.y","Day", "Group", "Group", "Cells.cm2.x6")
zoox.final$Day <- as.factor(zoox.final$Day)
zoox.final$Group <- as.factor(zoox.final$Group)
names(zoox.final)[1] <- "Fragment.ID"
write.csv(zoox.final, "../output/Physiology/Zoox.Calc.csv")

Chlorophyll a and c2

# Import data 
chla.raw.0708 <- read.csv("../data/Physiology/Chlorophyll/20210708_Chl.csv")
chla.raw.0716 <- read.csv("../data/Physiology/Chlorophyll/20210716_Chl.csv")

chla.meta <- read.csv("../data/Physiology/Chlorophyll/Chlorophyll_Meta.csv")


# Adding Run number to datasets
chla.raw.0708$Date <- 20210708
chla.raw.0716$Date <- 20210716

# Make a unique column 
chla.raw.0708$date.well <- paste(chla.raw.0708$Date, chla.raw.0708$Well, sep = "-")
chla.raw.0716$date.well <- paste(chla.raw.0716$Date, chla.raw.0716$Well, sep = "-")

chla.meta$date.well <- paste(chla.meta$Date, chla.meta$Well, sep = "-")

# Attaching vial metadata
chla.data.0708 <- merge(chla.raw.0708, chla.meta, by = "date.well")
chla.data.0716 <- merge(chla.raw.0716, chla.meta, by = "date.well")

# Blank 750nm correction for each run separately
Blank.0708 <- chla.data.0708 %>% 
  filter(Sample.Type == "Blank") %>%
  summarise(blk.avg = mean(Chl.750))

Blank.0716 <- chla.data.0716 %>% 
  filter(Sample.Type == "Blank") %>%
  summarise(blk.avg = mean(Chl.750))

# Subtracting 750 (blank) from 630 and 633 values, and accounting for path length (0.584 cm)
chla.data.0708$abs.630.corr <- (chla.data.0708$`Chl.630` - Blank.0708$blk.avg) / 0.584
chla.data.0708$abs.663.corr <- (chla.data.0708$`Chl.663` - Blank.0708$`blk.avg`) / 0.584

chla.data.0716$abs.630.corr <- (chla.data.0716$`Chl.630` - Blank.0716$blk.avg) / 0.584
chla.data.0716$abs.663.corr <- (chla.data.0716$`Chl.663` - Blank.0716$`blk.avg`) / 0.584

# Combining Datasets
chla.data.all <- rbind(chla.data.0708, chla.data.0716)

# Chlorp0hyll A concentration equation 
chla.data.all$chlA.ug.sample <- 11.43*chla.data.all$abs.663.corr - 0.64*chla.data.all$abs.630.corr

# Chlorophyll C2 concentration equation 
chla.data.all$chlC2.ug.sample <- 27.09*chla.data.all$abs.630.corr - 3.63*chla.data.all$abs.663.corr

# Attaching colony metadata
Chla.data.meta <- merge(chla.data.all, master.meta, by = "Vial")

# Standardization
Chla.data.meta$ChlA.ugcm2 <- (Chla.data.meta$chlA.ug.sample * Chla.data.meta$Homogenate_Vol.mL)/Chla.data.meta$Surface_Area_cm2 #Calculating concentration

Chla.data.meta$ChlC2.ugcm2 <- (Chla.data.meta$chlC2.ug.sample * Chla.data.meta$Homogenate_Vol.mL)/Chla.data.meta$Surface_Area_cm2 #Calculating concentration

# Removing well A9-20210716 because of pipette error
Chla.data.meta.clean <- Chla.data.meta %>%
  filter(date.well != "20210716-A9")

# Summarize per vial 
Chla.sum <- summarySE(Chla.data.meta.clean, measurevar="ChlA.ugcm2", groupvars=c("Vial", "Fragment.ID.y", "Day", "Group"))
Chla.sum2 <- Chla.sum %>%
  dplyr::select(Fragment.ID.y, Day, Group, ChlA.ugcm2)

ChlC2.sum <- summarySE(Chla.data.meta.clean, measurevar="ChlC2.ugcm2", groupvars=c("Vial", "Fragment.ID.y", "Day", "Group"))
ChlC2.sum2 <- ChlC2.sum %>%
  dplyr::select(Fragment.ID.y, Day, Group, ChlC2.ugcm2)

Chl.final <- merge(Chla.sum2, ChlC2.sum2, by = c("Fragment.ID.y", "Day", "Group"))

Chl.final$Day <- as.factor(Chl.final$Day)
Chl.final$Group <- as.factor(Chl.final$Group)
names(Chl.final)[1] <- "Fragment.ID"
write.csv(Chl.final, "../output/Physiology/Chl.Calc.csv")

Total Protein

TP.well.meta <- read.csv("../data/Physiology/Protein/Protein_Well_Meta.csv")

raw_20210701 <- read.csv("../data/Physiology/Protein/20210701_Protein.csv")
raw_20210702 <- read.csv("../data/Physiology/Protein/20210702_Protein.csv")
raw_20210706 <- read.csv("../data/Physiology/Protein/20210706_Protein.csv")
raw_20210707 <- read.csv("../data/Physiology/Protein/20210707_Protein.csv")

# Adding Run numbers into datasets
raw_20210701$run <- "20210701"
raw_20210702$run <- "20210702"
raw_20210706$run <- "20210706"
raw_20210707$run <- "20210707"

# Make a unique column for merging
raw_20210701$run.well <- paste(raw_20210701$run, raw_20210701$Well, sep = "-")
raw_20210702$run.well <- paste(raw_20210702$run, raw_20210702$Well, sep = "-")
raw_20210706$run.well <- paste(raw_20210706$run, raw_20210706$Well, sep = "-")
raw_20210707$run.well <- paste(raw_20210707$run, raw_20210707$Well, sep = "-")

TP.well.meta$run.well <- paste(TP.well.meta$Date, TP.well.meta$Well, sep = "-")

# Merge with metadata
TP.20210701 <- merge(raw_20210701, TP.well.meta, by = "run.well")
TP.20210702 <- merge(raw_20210702, TP.well.meta, by = "run.well")
TP.20210706 <- merge(raw_20210706, TP.well.meta, by = "run.well")
TP.20210707 <- merge(raw_20210707, TP.well.meta, by = "run.well")

# Subtract blanks means for each run 
TP.standardblank.01 <- TP.20210701 %>% 
  filter(Sample.Type == "Blank") %>%
  summarise(blk.avg = mean(X562))
TP.20210701$abs.corr <- TP.20210701$X562 - TP.standardblank.01$blk.avg

TP.standardblank.02 <- TP.20210702 %>% 
  filter(Sample.Type == "Blank") %>%
  summarise(blk.avg = mean(X562))
TP.20210702$abs.corr <- TP.20210702$X562 - TP.standardblank.02$blk.avg

TP.standardblank.06 <- TP.20210706 %>% 
  filter(Sample.Type == "Blank") %>%
  summarise(blk.avg = mean(X562))
TP.20210706$abs.corr <- TP.20210706$X562 - TP.standardblank.06$blk.avg

TP.standardblank.07 <- TP.20210707 %>% 
  filter(Sample.Type == "Blank") %>%
  summarise(blk.avg = mean(X562))
TP.20210707$abs.corr <- TP.20210707$X562 - TP.standardblank.07$blk.avg

# Run standards
TP.standard.01 <- TP.20210701 %>% 
  filter(Sample.Type == "Standard") 

TP.plot.S1<- ggplot(data = TP.standard.01, aes(x=Concentration, y=abs.corr))+
  ylab("Absorbance (nm)")+ xlab("Concentration") + 
  geom_point()+
  geom_smooth(method = "lm") +
  stat_regline_equation(label.y = 2.0, aes(label = ..eq.label..)) +
  stat_regline_equation(label.y = 1.75, aes(label = ..rr.label..)) +
  theme_bw() + theme(panel.border = element_blank(), panel.grid.major = element_blank(),
                     panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
TP.plot.S1

TP.lmstandard.01 <- lm (Concentration ~ abs.corr, data = TP.standard.01)
TP.lmsummary.01 <- summary(TP.lmstandard.01)


TP.standard.02 <- TP.20210702 %>% 
  filter(Sample.Type == "Standard") 

TP.plot.S2<- ggplot(data = TP.standard.02, aes(x=Concentration, y=abs.corr))+
  ylab("Absorbance (nm)")+ xlab("Concentration") + 
  geom_point()+
  geom_smooth(method = "lm") +
  stat_regline_equation(label.y = 2.0, aes(label = ..eq.label..)) +
  stat_regline_equation(label.y = 1.75, aes(label = ..rr.label..)) +
  theme_bw() + theme(panel.border = element_blank(), panel.grid.major = element_blank(),
                     panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
TP.plot.S2

TP.lmstandard.02 <- lm (Concentration ~ abs.corr, data = TP.standard.02)
TP.lmsummary.02 <- summary(TP.lmstandard.02)

TP.standard.06 <- TP.20210706 %>% 
  filter(Sample.Type == "Standard") 

TP.plot.S6<- ggplot(data = TP.standard.06, aes(x=Concentration, y=abs.corr))+
  ylab("Absorbance (nm)")+ xlab("Concentration") + 
  geom_point()+
  geom_smooth(method = "lm") +
  stat_regline_equation(label.y = 2.0, aes(label = ..eq.label..)) +
  stat_regline_equation(label.y = 1.75, aes(label = ..rr.label..)) +
  theme_bw() + theme(panel.border = element_blank(), panel.grid.major = element_blank(),
                     panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
TP.plot.S6

TP.lmstandard.06 <- lm (Concentration ~ abs.corr, data = TP.standard.06)
TP.lmsummary.06 <- summary(TP.lmstandard.06)


TP.standard.07 <- TP.20210707 %>% 
  filter(Sample.Type == "Standard") 

TP.plot.S7<- ggplot(data = TP.standard.07, aes(x=Concentration, y=abs.corr))+
  ylab("Absorbance (nm)")+ xlab("Concentration") + 
  geom_point()+
  geom_smooth(method = "lm") +
  stat_regline_equation(label.y = 2.0, aes(label = ..eq.label..)) +
  stat_regline_equation(label.y = 1.75, aes(label = ..rr.label..)) +
  theme_bw() + theme(panel.border = element_blank(), panel.grid.major = element_blank(),
                     panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
TP.plot.S7

TP.lmstandard.07 <- lm (Concentration ~ abs.corr, data = TP.standard.07)
TP.lmsummary.07 <- summary(TP.lmstandard.07)


# Extrapolate concentration values

TP.Sample.01 <- TP.20210701 %>% #subsetting Samples
  filter(Sample.Type == "Sample") 
TP.Sample.01$ConcentrationS <- predict(TP.lmstandard.01, newdata = TP.Sample.01) #using model to get concentration

TP.Sample.02 <- TP.20210702 %>% #subsetting Samples
  filter(Sample.Type == "Sample") 
TP.Sample.02$ConcentrationS <- predict(TP.lmstandard.02, newdata = TP.Sample.02) #using model to get concentration

TP.Sample.06 <- TP.20210706 %>% #subsetting Samples
  filter(Sample.Type == "Sample") 
TP.Sample.06$ConcentrationS <- predict(TP.lmstandard.06, newdata = TP.Sample.06) #using model to get concentration

TP.Sample.07 <- TP.20210707 %>% #subsetting Samples
  filter(Sample.Type == "Sample") 
TP.Sample.07$ConcentrationS <- predict(TP.lmstandard.07, newdata = TP.Sample.07) #using model to get concentration

# Row Bind DFs together
TP.all <- rbind(TP.Sample.01, TP.Sample.02, TP.Sample.06, TP.Sample.07)

# Adding metadata
TP.all.meta <- merge(TP.all, master.meta, by = "Vial")

# Standardization 
TP.all.meta$Total.Protein.ugcm2 <- (TP.all.meta$ConcentrationS * TP.all.meta$Homogenate_Vol.mL*1.45)/TP.all.meta$Surface_Area_cm2 #Calculating concentration. 1.45 = Dilution factor of acid+base in sample
TP.all.meta$Total.Protein.mgcm2 <- TP.all.meta$Total.Protein.ugcm2/1000 

# Summarize per vial 
TP.sum <- summarySE(TP.all.meta, measurevar="Total.Protein.mgcm2", groupvars=c("Vial", "Fragment.ID.y", "Day", "Group", "Type"))
TP.sum2 <- TP.sum %>%
  dplyr::select(Fragment.ID.y, Day, Group, Type, Total.Protein.mgcm2)

# Converting long to wide format
TP.final<- TP.sum2 %>% 
  pivot_wider(names_from = Type, values_from = Total.Protein.mgcm2)

# Calculating Host to Symbiont Ratio
TP.final$HostSymRatioProtein <- TP.final$Coral / TP.final$Symbiont

# Renaming columns 

TP.final2 <- TP.final %>%
  dplyr::rename(Coral_Protein_mgcm2 = Coral, 
         Symbiont_Protein_mgcm2 = Symbiont)

TP.final2$Day <- as.factor(TP.final2$Day)
TP.final2$Group <- as.factor(TP.final2$Group)
names(TP.final2)[1] <- "Fragment.ID"
write.csv(TP.final2, "../output/Physiology/Protein.Calc.csv")

Carbohydrate

# Read in datafiles

carb.data.0715 <- read.csv("../data/Physiology/Carbohydrate/20210715_Carb.csv")
carb.data.0716 <- read.csv("../data/Physiology/Carbohydrate/20210716_Carb.csv")
carb.data.0719_1 <- read.csv("../data/Physiology/Carbohydrate/20210719-1_Carb.csv")
carb.data.0719_2 <- read.csv("../data/Physiology/Carbohydrate/20210719-2_Carb.csv")

carb.meta <- read.csv("../data/Physiology/Carbohydrate/Carb_Meta.csv")

# Adding a date column

carb.data.0715$Date <- "20210715"
carb.data.0716$Date <- "20210716"
carb.data.0719_1$Date <- "20210719-1"
carb.data.0719_2$Date <- "20210719-2"
  
# Making unique columns by date/well

carb.data.0715$date.well <- paste(carb.data.0715$Date, carb.data.0715$Well, sep = "-")
carb.data.0716$date.well <- paste(carb.data.0716$Date, carb.data.0716$Well, sep = "-")
carb.data.0719_1$date.well <- paste(carb.data.0719_1$Date, carb.data.0719_1$Well, sep = "-")
carb.data.0719_2$date.well <- paste(carb.data.0719_2$Date, carb.data.0719_2$Well, sep = "-")

carb.meta$date.well <- paste(carb.meta$Date, carb.meta$Well, sep = "-")

# Merging Metadata 

carb.data.meta.0715 <- merge(carb.data.0715, carb.meta, by = "date.well") #has standard
carb.data.meta.0716 <- merge(carb.data.0716, carb.meta, by = "date.well") #use 0715 standard
carb.data.meta.0719_1 <- merge(carb.data.0719_1, carb.meta, by = "date.well") #has standard
carb.data.meta.0719_2 <- merge(carb.data.0719_2, carb.meta, by = "date.well") #use 0719-1 standard

# Blank correction for each corresponding run 

Blank.0715 <- carb.data.meta.0715 %>% 
  filter(Sample.Type == "Blank") %>%
  summarise(blk.avg = mean(X485))

carb.data.meta.0715$abs.corr <- carb.data.meta.0715$X485 - Blank.0715$blk.avg
carb.data.meta.0716$abs.corr <- carb.data.meta.0716$X485 - Blank.0715$blk.avg

Blank.0719 <- carb.data.meta.0719_1 %>% 
  filter(Sample.Type == "Blank") %>%
  summarise(blk.avg = mean(X485))

carb.data.meta.0719_1$abs.corr <- carb.data.meta.0719_1$X485 - Blank.0719$blk.avg
carb.data.meta.0719_2$abs.corr <- carb.data.meta.0719_2$X485 - Blank.0719$blk.avg


# Standard curve 

Standard.0715 <- carb.data.meta.0715 %>% 
  filter(Sample.Type == "Standard") 

Standard.0715.plot <- ggplot(data = Standard.0715, aes(x=Concentration, y=abs.corr))+
  ylab("Absorbance (nm)")+ xlab("Carbohydrate (mg/mL)") + 
  geom_point()+
  geom_smooth(method = "lm") +
  stat_regline_equation(label.y = 1.0, aes(label = ..eq.label..)) +
  stat_regline_equation(label.y = 0.75, aes(label = ..rr.label..)) +
  theme_bw() + theme(panel.border = element_blank(), panel.grid.major = element_blank(),
                     panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
Standard.0715.plot

lmstandard.0715 <- lm (Concentration ~ abs.corr, data = Standard.0715)
lmsummary.0715 <- summary(lmstandard.0715) 


Standard.0719 <- carb.data.meta.0719_1 %>% 
  filter(Sample.Type == "Standard") 

Standard.0719.plot <- ggplot(data = Standard.0719, aes(x=Concentration, y=abs.corr))+
  ylab("Absorbance (nm)")+ xlab("Carbohydrate (mg/mL)") + 
  geom_point()+
  geom_smooth(method = "lm") +
  stat_regline_equation(label.y = 1.0, aes(label = ..eq.label..)) +
  stat_regline_equation(label.y = 0.75, aes(label = ..rr.label..)) +
  theme_bw() + theme(panel.border = element_blank(), panel.grid.major = element_blank(),
                     panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
Standard.0719.plot

lmstandard.0719 <- lm (Concentration ~ abs.corr, data = Standard.0719)
lmsummary.0719 <- summary(lmstandard.0719) 

# Obtaining concentration values from standard curve

Samples.0715 <- carb.data.meta.0715 %>% #subsetting Samples
  filter(Sample.Type == "Sample") 
Samples.0715$Concentration <- predict(lmstandard.0715, newdata = Samples.0715) #using model to get concentration

Samples.0716 <- carb.data.meta.0716 %>% #subsetting Samples
  filter(Sample.Type == "Sample") 
Samples.0716$Concentration <- predict(lmstandard.0715, newdata = Samples.0716) #using model to get concentration

Samples.0719_1 <- carb.data.meta.0719_1 %>% #subsetting Samples
  filter(Sample.Type == "Sample") 
Samples.0719_1$Concentration <- predict(lmstandard.0719, newdata = Samples.0719_1) #using model to get concentration

Samples.0719_2 <- carb.data.meta.0719_2 %>% #subsetting Samples
  filter(Sample.Type == "Sample") 
Samples.0719_2$Concentration <- predict(lmstandard.0719, newdata = Samples.0719_2) #using model to get concentration

# Merging metadata 

Samples.carb.all <- rbind(Samples.0715, Samples.0716, Samples.0719_1, Samples.0719_2)

samples.carb.meta <- merge(Samples.carb.all, master.meta, by = "Vial")

# Accounting for dilution factor (1000/10) and Normalizing to homogenate volume and surface area

samples.carb.meta$Carb.mgcm2 <- (samples.carb.meta$Concentration * samples.carb.meta$Homogenate_Vol.mL * 10)/samples.carb.meta$Surface_Area_cm2 

#Summarize per vial 
carb.sum <- summarySE(samples.carb.meta, measurevar="Carb.mgcm2", groupvars=c("Vial", "Fragment.ID.y", "Day", "Group", "Type"))
carb.sum2 <- carb.sum %>%
  dplyr::select(Fragment.ID.y, Day, Group, Type, Carb.mgcm2)

# Converting long to wide format
Carb.final<- carb.sum2 %>% 
  pivot_wider(names_from = Type, values_from = Carb.mgcm2)

# Calculating Host to Symbiont Ratio
Carb.final$HostSymRatioCarb <- Carb.final$Coral / Carb.final$Symbiont

# Renaming columns 

Carb.final2 <- Carb.final %>%
  dplyr::rename(Coral_Carb_mgcm2 = Coral, 
                Symbiont_Carb_mgcm2 = Symbiont)

Carb.final2$Day <- as.factor(Carb.final2$Day)
Carb.final2$Group <- as.factor(Carb.final2$Group)
names(Carb.final2)[1] <- "Fragment.ID"
write.csv(Carb.final2, "../output/Physiology/Carb.Calc.csv")

Master Datasheet and Normalizations

resp_dataset <- read.csv("../output/Physiology/Resp.Calc.csv")
carb_dataset <- read.csv("../output/Physiology/Carb.Calc.csv")
prot_dataset <- read.csv("../output/Physiology/Protein.Calc.csv")
sym_dataset <- read.csv("../output/Physiology/Zoox.Calc.csv")
chl_dataset <- read.csv("../output/Physiology/Chl.Calc.csv")
frag.meta <- read.csv("../data/Physiology/Airbrush_Metadata.csv")

#selecting SA column 

frag.sa <- frag.meta %>%
  dplyr::select(Fragment.ID, Day, Group, Surface_Area_cm2)

# Joining datasets

dfs <- list(resp_dataset,
            carb_dataset,
            prot_dataset,
            sym_dataset,
            chl_dataset, 
            frag.sa)

master <- join_all(dfs, by=c("Fragment.ID", "Day", "Group")) %>%
  dplyr::select(-ends_with("X"))

# Normalization to per sym cell 

master$Cells <- master$Cells.cm2.x6 * master$Surface_Area_cm2 * 1000000 # getting absolute number of cells 
master$ChlA.ugcell <- ((master$ChlA.ugcm2*master$Surface_Area_cm2) / (master$Cells)) 
master$ChlC2.ugcell <- ((master$ChlC2.ugcm2*master$Surface_Area_cm2) / (master$Cells))
master$Carb.mgcell <- ((master$Symbiont_Carb_mgcm2*master$Surface_Area_cm2) / (master$Cells))
master$Protein.mgcell <- ((master$Symbiont_Protein_mgcm2*master$Surface_Area_cm2) / (master$Cells))

  
master$Day <- as.factor(master$Day)
master$Group <- as.factor(master$Group)
master$Fragment.ID <- as.factor(master$Fragment.ID)

master
##    Fragment.ID Day     Group Pnet_umol.cm2.hr Pgross_umol.cm2.hr
## 1          R11  52 Mortality      -0.23189241       0.1020934734
## 2          R11   0 Mortality       0.50582731       1.1868464067
## 3          R11  37 Mortality      -0.36245920       0.0033265282
## 4          R17  52   Control       0.26802023       0.6920767697
## 5          R17  37   Control       0.71085019       1.3291667409
## 6          R17   0   Control       0.41877850       1.1186333864
## 7          R19  52  Bleached      -0.41907750       0.0567272258
## 8          R19   0  Bleached       0.38354936       1.1095961967
## 9          R19  37  Bleached      -0.11848265       0.1411971380
## 10         R20  52  Bleached      -0.35947860      -0.0003639834
## 11         R20  37  Bleached      -0.16116833       0.1184719656
## 12         R20   0  Bleached       0.16410711       0.6424367190
## 13         R23  52   Control      -0.04329761       0.3978362514
## 14         R23   0   Control       0.24238553       0.7307221520
## 15         R23  37   Control       0.23486158       0.8538006108
## 16         R26  52 Mortality      -0.16624709       0.2905505959
## 17         R26  37 Mortality      -0.32821468       0.1267031397
## 18         R26   0 Mortality       0.15241072       0.7662224634
## 19         R28  52 Mortality      -0.19605011       0.3981000459
## 20         R28   0 Mortality       0.21938199       1.0187704774
## 21         R28  37 Mortality      -0.19076104       0.1604887254
## 22         R29  52  Bleached      -0.20851651       0.1901868201
## 23         R29  37  Bleached      -0.18425232       0.1143434084
## 24         R29   0  Bleached       0.32886894       1.0688889329
## 25         R32  37   Control       0.19991505       0.5733927695
## 26         R32  52   Control       0.12338534       0.5100397489
## 27         R32   0   Control       0.03985412       0.7318985087
## 28         R35  52 Mortality      -0.24541776       0.2131728043
## 29         R35  37 Mortality      -0.17148385       0.2335203154
## 30         R35   0 Mortality       0.02967193       0.8438741572
## 31         R36  52 Mortality      -0.37237629       0.0934852282
## 32         R36  37 Mortality      -0.44603538      -0.0156581586
## 33         R36   0 Mortality       0.08499423       0.6582026486
## 34         R37  37  Bleached      -0.34553509       0.1061178335
## 35         R37  52  Bleached      -0.30877729       0.0663087355
## 36         R37   0  Bleached       0.25089923       1.1171429486
## 37         R40  37   Control       0.30772274       0.6987146249
## 38         R40   0   Control       0.22809026       0.9602083084
## 39         R40  52   Control       0.01742511       0.5202652851
## 40          R7  37   Control       0.73880159       1.2836098810
## 41          R7   0   Control       0.15738175       0.8461748059
## 42          R7  52   Control       0.15520373       0.7734739299
## 43          R8  37  Bleached      -0.44205832       0.1024291844
## 44          R8   0  Bleached       0.26755533       0.7632090733
## 45          R8  52  Bleached      -0.28186064       0.0627559505
##    abs.Rdark_umol.cm2.hr           PR Coral_Carb_mgcm2 Symbiont_Carb_mgcm2
## 1              0.3339859  0.305682003         6.849361           1.7643484
## 2              0.6810191  1.742750543        11.308695           3.1202253
## 3              0.3657857  0.009094199        12.095280           3.8631560
## 4              0.4240565  1.632038911         6.839514           2.9490419
## 5              0.6183166  2.149654148         5.563662           3.0202852
## 6              0.6998549  1.598379057         7.566091           2.7553698
## 7              0.4758047  0.119223754        10.949188           3.7515966
## 8              0.7260468  1.528270817         9.239969           2.1188079
## 9              0.2596798  0.543735581        11.808212           3.0126430
## 10             0.3591146 -0.001013558         3.501538           2.4675249
## 11             0.2796403  0.423658414        11.267394           3.7730921
## 12             0.4783296  1.343083741        13.414633           3.3379630
## 13             0.4411339  0.901849265         7.433731           3.1970812
## 14             0.4883366  1.496349292         7.162866           3.0593058
## 15             0.6189390  1.379458335         9.196945           3.1675563
## 16             0.4567977  0.636059691         4.323523           2.5349244
## 17             0.4549178  0.278518741         9.142666           3.8381270
## 18             0.6138117  1.248302056         8.903752           3.2006537
## 19             0.5941502  0.670032725         5.728313           2.3687238
## 20             0.7993885  1.274437258         9.784738           1.3047815
## 21             0.3512498  0.456907706         6.014000           2.1604203
## 22             0.3987033  0.477013371         8.957309           2.5973352
## 23             0.2985957  0.382937186         8.994702           2.0654242
## 24             0.7400200  1.444405471        11.828278           3.2892179
## 25             0.3734777  1.535279733        11.895859           4.2947185
## 26             0.3866544  1.319110125         7.044523           2.7746750
## 27             0.6920444  1.057588969        11.424910           4.5526846
## 28             0.4585906  0.464843411         7.295854           3.8457757
## 29             0.4050042  0.576587447         5.638928           3.6298139
## 30             0.8142022  1.036442949        12.973316           7.4011411
## 31             0.4658615  0.200671712         7.707399           3.0474493
## 32             0.4303772 -0.036382406         7.803940           0.4781953
## 33             0.5732084  1.148278058        10.560248           2.6024627
## 34             0.4516529  0.234954384        13.123489           5.3589500
## 35             0.3750860  0.176782740         6.351416           2.0801262
## 36             0.8662437  1.289640462        10.656405           3.3870259
## 37             0.3909919  1.787031026         7.585425           2.6984209
## 38             0.7321180  1.311548472         7.970075           3.6814867
## 39             0.5028402  1.034653372         7.523901           3.3213336
## 40             0.5448083  2.356076265         9.237811           4.3413046
## 41             0.6887931  1.228489171         6.444712           2.4862677
## 42             0.6182702  1.251028970         6.673854           2.6173143
## 43             0.5444875  0.188120358         5.671352           2.4961686
## 44             0.4956537  1.539802914        10.097315           3.7617888
## 45             0.3446166  0.182103684         6.963212           3.5262658
##    HostSymRatioCarb Coral_Protein_mgcm2 Symbiont_Protein_mgcm2
## 1          3.882091           0.7777801              1.0472301
## 2          3.624320           1.2378482              1.7432650
## 3          3.130932           1.8390078              1.5447205
## 4          2.319233           0.8736717              1.1589387
## 5          1.842098           0.6569932              1.0931865
## 6          2.745944           1.1305265              1.2609914
## 7          2.918541           0.8252890              1.2452956
## 8          4.360928           0.9373076              0.8384386
## 9          3.919552           1.1680675              0.9544774
## 10         1.419049           1.1041223              1.4536745
## 11         2.986249           1.5541780              1.6921869
## 12         4.018808           1.4416994              1.4976326
## 13         2.325162           0.8660496              1.1289232
## 14         2.341337           1.0767521              1.2553001
## 15         2.903483           0.9825597              1.3234874
## 16         1.705583           0.7562860              1.2217228
## 17         2.382065           1.5382525              1.9977511
## 18         2.781854           1.0429995              1.4323859
## 19         2.418312           0.5977390              0.8906722
## 20         7.499139           1.1932539              1.3533129
## 21         2.783718           0.9747326              0.9724740
## 22         3.448653           1.1988134              1.5509649
## 23         4.354893           1.1360409              1.2085604
## 24         3.596076           1.5561040              1.6477717
## 25         2.769881           1.4273253              1.4083638
## 26         2.538864           1.1112201              1.3982868
## 27         2.509489           1.7837592              2.2583977
## 28         1.897109           1.2431900              1.7403353
## 29         1.553503           1.1883370              1.3598986
## 30         1.752881           1.6339445              1.3969345
## 31         2.529131           1.1020099              1.4598403
## 32        16.319564           0.9821318              0.8657490
## 33         4.057790           1.1750891              1.1892187
## 34         2.448892           1.6609407              2.0817057
## 35         3.053380           1.1211654              1.2873024
## 36         3.146242           1.7877324              1.7211346
## 37         2.811061           0.8884828              1.1446421
## 38         2.164907           0.9103872              1.0356749
## 39         2.265325           0.9648052              1.4423972
## 40         2.127888           0.8450959              1.1170357
## 41         2.592123           0.9654668              1.3912152
## 42         2.549886           0.6772846              1.0224350
## 43         2.272023           0.8731718              1.2980929
## 44         2.684179           1.0600605              1.5346158
## 45         1.974670           0.8337565              1.2834567
##    HostSymRatioProtein Cells.cm2.x6 ChlA.ugcm2 ChlC2.ugcm2 Surface_Area_cm2
## 1            0.7427022    0.1540965  0.3275101  0.09132213           10.999
## 2            0.7100746    1.4437051  7.4072233  2.48951700            4.358
## 3            1.1905117    0.4190049  1.3837000  0.49063696            2.970
## 4            0.7538550    0.8260314  2.6574699  0.83689365            5.478
## 5            0.6009892    1.3826249  5.6995893  1.64190686            8.869
## 6            0.8965378    2.0221236 10.3365517  3.33457116            9.462
## 7            0.6627254    0.1963828  1.3405398  0.38270509            4.366
## 8            1.1179204    1.1433489  4.2424029  1.36950514            8.134
## 9            1.2237769    0.3792668  1.0811195  1.01473339            2.373
## 10           0.7595389    0.2994328  0.9156895  1.76202006            7.457
## 11           0.9184435    0.4793840  3.5937723  1.05426883            3.355
## 12           0.9626523    1.6895325  8.2288377  2.63920650            4.346
## 13           0.7671466    0.7009038  4.8898627  5.37157161            5.053
## 14           0.8577646    1.3687119  5.1644179  1.55973806            5.851
## 15           0.7424020    1.1838762  2.0581962  0.25032277            4.329
## 16           0.6190324    0.4080261  1.2990092  0.27255717           10.416
## 17           0.7699920    0.4046492  1.5966270  0.50438075            4.438
## 18           0.7281554    1.3479952  6.0547269  2.04455744            5.570
## 19           0.6711100    0.1033754  0.7077547  0.17061524            8.733
## 20           0.8817280    1.6370114  9.2615161  7.27550581            8.374
## 21           1.0023225    0.1921334  0.5019107  0.21612215            6.130
## 22           0.7729468    0.2255075  0.2815202  0.02862173            5.851
## 23           0.9399951    0.1618713  0.9437220  0.93641702            5.028
## 24           0.9443687    1.0647622  3.5832597  1.41391240            3.518
## 25           1.0134635    1.0270481  4.1029036  1.38907203            2.779
## 26           0.7947011    1.1101703  2.7164310  1.05396663            9.448
## 27           0.7898340    2.4036319 13.8831786  4.10523954            2.244
## 28           0.7143393    0.8808811  1.2484700  0.60104648            4.846
## 29           0.8738423    0.3630730  1.3484589  0.46397694            3.870
## 30           1.1696644    1.8293839  6.1702672  1.85880997            2.110
## 31           0.7548839    0.3418154  0.7865949  0.27474975            5.266
## 32           1.1344302    0.4723115  0.7238693  0.27262691            5.497
## 33           0.9881186    1.3758214  4.9521489  1.59004149            3.348
## 34           0.7978749    0.5588378  4.1108403  1.12810839            2.144
## 35           0.8709418    0.1145972  0.5022336  0.14783410            7.490
## 36           1.0386941    1.4702543  5.2101219  1.63831837            3.670
## 37           0.7762101    0.8459411  3.5594404  1.15392072            8.196
## 38           0.8790280    1.7612390  9.0560350  2.84515065            1.824
## 39           0.6688900    0.9430216  2.6613505  0.80002162            5.037
## 40           0.7565523    1.1229514  5.1353209  1.67873098            4.434
## 41           0.6939737    1.2677780  4.9487060  1.40710187            8.578
## 42           0.6624232    0.9563789  2.0493608  0.73343085           10.809
## 43           0.6726574    0.1335096  1.0028003  0.39763902            6.429
## 44           0.6907660    2.1099318 11.5075025  3.63966661            4.334
## 45           0.6496180    0.2329879  1.3589985  0.44224326            5.532
##         Cells  ChlA.ugcell ChlC2.ugcell  Carb.mgcell Protein.mgcell
## 1   1694907.4 2.125357e-06 5.926295e-07 1.144963e-05   6.795937e-06
## 2   6291666.7 5.130704e-06 1.724394e-06 2.161262e-06   1.207494e-06
## 3   1244444.4 3.302348e-06 1.170958e-06 9.219836e-06   3.686641e-06
## 4   4525000.0 3.217154e-06 1.013150e-06 3.570133e-06   1.403020e-06
## 5  12262500.0 4.122296e-06 1.187529e-06 2.184457e-06   7.906602e-07
## 6  19133333.3 5.111731e-06 1.649044e-06 1.362612e-06   6.235976e-07
## 7    857407.4 6.826156e-06 1.948771e-06 1.910349e-05   6.341164e-06
## 8   9300000.0 3.710506e-06 1.197802e-06 1.853159e-06   7.333182e-07
## 9    900000.0 2.850552e-06 2.675514e-06 7.943335e-06   2.516639e-06
## 10  2232870.4 3.058080e-06 5.884526e-06 8.240663e-06   4.854760e-06
## 11  1608333.3 7.496646e-06 2.199216e-06 7.870709e-06   3.529919e-06
## 12  7342708.3 4.870482e-06 1.562093e-06 1.975673e-06   8.864183e-07
## 13  3541666.7 6.976511e-06 7.663779e-06 4.561370e-06   1.610668e-06
## 14  8008333.3 3.773196e-06 1.139566e-06 2.235171e-06   9.171398e-07
## 15  5125000.0 1.738523e-06 2.114434e-07 2.675581e-06   1.117927e-06
## 16  4250000.0 3.183642e-06 6.679895e-07 6.212652e-06   2.994227e-06
## 17  1795833.3 3.945706e-06 1.246464e-06 9.485072e-06   4.936995e-06
## 18  7508333.3 4.491653e-06 1.516739e-06 2.374381e-06   1.062605e-06
## 19   902777.8 6.846449e-06 1.650443e-06 2.291379e-05   8.615897e-06
## 20 13708333.3 5.657576e-06 4.444383e-06 7.970510e-07   8.266973e-07
## 21  1177777.8 2.612303e-06 1.124855e-06 1.124438e-05   5.061452e-06
## 22  1319444.4 1.248385e-06 1.269214e-07 1.151773e-05   6.877664e-06
## 23   813888.9 5.830076e-06 5.784948e-06 1.275967e-05   7.466181e-06
## 24  3745833.3 3.365315e-06 1.327914e-06 3.089157e-06   1.547549e-06
## 25  2854166.7 3.994851e-06 1.352490e-06 4.181614e-06   1.371273e-06
## 26 10488888.9 2.446860e-06 9.493738e-07 2.499324e-06   1.259525e-06
## 27  5393750.0 5.775917e-06 1.707932e-06 1.894086e-06   9.395772e-07
## 28  4268750.0 1.417297e-06 6.823242e-07 4.365828e-06   1.975676e-06
## 29  1405092.6 3.714016e-06 1.277916e-06 9.997476e-06   3.745524e-06
## 30  3860000.0 3.372866e-06 1.016085e-06 4.045701e-06   7.636093e-07
## 31  1800000.0 2.301227e-06 8.037957e-07 8.915482e-06   4.270844e-06
## 32  2596296.3 1.532610e-06 5.772185e-07 1.012458e-06   1.833004e-06
## 33  4606250.0 3.599413e-06 1.155703e-06 1.891570e-06   8.643699e-07
## 34  1198148.1 7.356053e-06 2.018669e-06 9.589456e-06   3.725063e-06
## 35   858333.3 4.382597e-06 1.290032e-06 1.815163e-05   1.123328e-05
## 36  5395833.3 3.543688e-06 1.114310e-06 2.303701e-06   1.170637e-06
## 37  6933333.3 4.207669e-06 1.364067e-06 3.189845e-06   1.353099e-06
## 38  3212500.0 5.141855e-06 1.615426e-06 2.090282e-06   5.880377e-07
## 39  4750000.0 2.822152e-06 8.483598e-07 3.522012e-06   1.529548e-06
## 40  4979166.7 4.573057e-06 1.494927e-06 3.865977e-06   9.947320e-07
## 41 10875000.0 3.903448e-06 1.109896e-06 1.961122e-06   1.097365e-06
## 42 10337500.0 2.142833e-06 7.668831e-07 2.736692e-06   1.069069e-06
## 43   858333.3 7.511072e-06 2.978355e-06 1.869655e-05   9.722842e-06
## 44  9144444.4 5.453969e-06 1.725016e-06 1.782896e-06   7.273296e-07
## 45  1288888.9 5.832915e-06 1.898139e-06 1.513498e-05   5.508685e-06

Boxplots and statistics

###############################

TP_Coral_Box <- ggplot(master, aes(x=Day, y=Coral_Protein_mgcm2, fill = Group)) +
  geom_boxplot(width=.5, outlier.shape= NA, position = position_dodge(width = 0.5)) +
  stat_summary(fun=mean, geom="line", aes(group=Group, color = Group), position = position_dodge(width = 0.5))  + 
  #  stat_summary(fun=mean, geom="point")
  geom_point(pch = 21, position=position_jitterdodge(dodge.width=0.5)) +
  #  ylim(0,0.5) +
  scale_fill_manual(values=c("#E7B800", "#00AFBB", "#FC4E07")) +
  scale_color_manual(values=c("#E7B800", "#00AFBB", "#FC4E07")) + 
  scale_x_discrete(labels=c("0" = "Day 0", "37" = "Day 37", "52" = "Day 52")) +
  xlab("Timepoint") +  ylab(expression("Total Protein " (mg~cm^{-2}))) + #Axis titles
  theme_bw() + theme(panel.border = element_rect(color="black", fill=NA, size=0.75), panel.grid.major = element_blank(), #Makes background theme white
                     panel.grid.minor = element_blank(), axis.line = element_blank()) +
  theme(axis.text = element_text(size = 15, color = "black"),
        axis.title = element_text(size = 18, color = "black")) + ggtitle("Coral Total Protein")
TP_Coral_Box

TP_Coral_lmer <- lmer(Coral_Protein_mgcm2~Group*Day+(1|Fragment.ID), data = master)
# qqnorm(resid(TP_Coral_lmer))
# qqline(resid(TP_Coral_lmer))
# 
# boxplot(resid(TP_Coral_lmer)~master$Group)
# boxplot(resid(TP_Coral_lmer)~master$Day)
anova(TP_Coral_lmer, type="II")
## Type II Analysis of Variance Table with Satterthwaite's method
##            Sum Sq Mean Sq NumDF DenDF F value   Pr(>F)   
## Group     0.10064 0.05032     2    12  1.0667 0.374642   
## Day       0.86013 0.43006     2    24  9.1168 0.001134 **
## Group:Day 0.16510 0.04127     4    24  0.8750 0.493418   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#capture.output(anova(TP_Coral_lmer, type="II"), file = "../output/Statistics/TP_Coral_lmer.csv")

###############################

TP_Sym_SA_Box <- ggplot(master, aes(x=Day, y=Symbiont_Protein_mgcm2, fill = Group)) +
  geom_boxplot(width=.5, outlier.shape= NA, position = position_dodge(width = 0.5)) +
  stat_summary(fun=mean, geom="line", aes(group=Group, color = Group), position = position_dodge(width = 0.5))  + 
  #  stat_summary(fun=mean, geom="point")
  geom_point(pch = 21, position=position_jitterdodge(dodge.width=0.5)) +
  #  ylim(0,0.5) +
  scale_fill_manual(values=c("#E7B800", "#00AFBB", "#FC4E07")) +
  scale_color_manual(values=c("#E7B800", "#00AFBB", "#FC4E07")) + 
  scale_x_discrete(labels=c("0" = "Day 0", "37" = "Day 37", "52" = "Day 52")) +
  xlab("Timepoint") +  ylab(expression("Total Protein " (mg~cm^{-2}))) + #Axis titles
  theme_bw() + theme(panel.border = element_rect(color="black", fill=NA, size=0.75), panel.grid.major = element_blank(), #Makes background theme white
                     panel.grid.minor = element_blank(), axis.line = element_blank()) +
  theme(axis.text = element_text(size = 15, color = "black"),
        axis.title = element_text(size = 18, color = "black")) + ggtitle("Symbiont Total Protein (per cm2)")
TP_Sym_SA_Box

TP_Sym_lmer <- lmer(Symbiont_Protein_mgcm2~Group*Day+(1|Fragment.ID), data = master)
# qqnorm(resid(TP_Sym_lmer))
# qqline(resid(TP_Sym_lmer))
# 
# boxplot(resid(TP_Sym_lmer)~master$Group)
# boxplot(resid(TP_Sym_lmer)~master$Day)
anova(TP_Sym_lmer, type="II")
## Type II Analysis of Variance Table with Satterthwaite's method
##             Sum Sq  Mean Sq NumDF DenDF F value Pr(>F)
## Group     0.060275 0.030138     2    12  0.3625 0.7033
## Day       0.171463 0.085732     2    24  1.0312 0.3718
## Group:Day 0.065453 0.016363     4    24  0.1968 0.9376
#capture.output(anova(TP_Sym_lmer, type="II"), file = "../output/Statistics/TP_Sym_lmer.csv")

###############################

TP_Sym_Cell_Box <- ggplot(master, aes(x=Day, y=Protein.mgcell, fill = Group)) +
  geom_boxplot(width=.5, outlier.shape= NA, position = position_dodge(width = 0.5)) +
  stat_summary(fun=mean, geom="line", aes(group=Group, color = Group), position = position_dodge(width = 0.5))  + 
  #  stat_summary(fun=mean, geom="point")
  geom_point(pch = 21, position=position_jitterdodge(dodge.width=0.5)) +
  #  ylim(0,0.5) +
  scale_fill_manual(values=c("#E7B800", "#00AFBB", "#FC4E07")) +
  scale_color_manual(values=c("#E7B800", "#00AFBB", "#FC4E07")) + 
  scale_x_discrete(labels=c("0" = "Day 0", "37" = "Day 37", "52" = "Day 52")) +
  xlab("Timepoint") +  ylab(expression("Total Protein " (mg~cell^{-1}))) + #Axis titles
  theme_bw() + theme(panel.border = element_rect(color="black", fill=NA, size=0.75), panel.grid.major = element_blank(), #Makes background theme white
                     panel.grid.minor = element_blank(), axis.line = element_blank()) +
  theme(axis.text = element_text(size = 15, color = "black"),
        axis.title = element_text(size = 18, color = "black")) + ggtitle("Symbiont Total Protein (per cell)")
TP_Sym_Cell_Box

TP_Sym_lmer <- lmer(Symbiont_Protein_mgcm2~Group*Day+(1|Fragment.ID), data = master)
# qqnorm(resid(TP_Sym_lmer))
# qqline(resid(TP_Sym_lmer))
# 
# boxplot(resid(TP_Sym_lmer)~master$Group)
# boxplot(resid(TP_Sym_lmer)~master$Day)
anova(TP_Sym_lmer, type="II")
## Type II Analysis of Variance Table with Satterthwaite's method
##             Sum Sq  Mean Sq NumDF DenDF F value Pr(>F)
## Group     0.060275 0.030138     2    12  0.3625 0.7033
## Day       0.171463 0.085732     2    24  1.0312 0.3718
## Group:Day 0.065453 0.016363     4    24  0.1968 0.9376
#capture.output(anova(TP_Sym_lmer, type="II"), file = "../output/Statistics/TP_Sym_lmer.csv")

###############################

TP_CS_Box <- ggplot(master, aes(x=Day, y=HostSymRatioProtein, fill = Group)) +
  geom_boxplot(width=.5, outlier.shape= NA, position = position_dodge(width = 0.5)) +
  stat_summary(fun=mean, geom="line", aes(group=Group, color = Group), position = position_dodge(width = 0.5))  + 
  #  stat_summary(fun=mean, geom="point")
  geom_point(pch = 21, position=position_jitterdodge(dodge.width=0.5)) +
  #  ylim(0,0.5) +
  scale_fill_manual(values=c("#E7B800", "#00AFBB", "#FC4E07")) +
  scale_color_manual(values=c("#E7B800", "#00AFBB", "#FC4E07")) + 
  scale_x_discrete(labels=c("0" = "Day 0", "37" = "Day 37", "52" = "Day 52")) +
  xlab("Timepoint") +  ylab("Coral to Symbiont Ratio (Total Protein)") + #Axis titles
  theme_bw() + theme(panel.border = element_rect(color="black", fill=NA, size=0.75), panel.grid.major = element_blank(), #Makes background theme white
                     panel.grid.minor = element_blank(), axis.line = element_blank()) +
  theme(axis.text = element_text(size = 15, color = "black"),
        axis.title = element_text(size = 18, color = "black")) + ggtitle("Host to Symbiont Ratio of Total Protein")
TP_CS_Box

TP_CS_lmer <- lmer(HostSymRatioProtein~Group*Day+(1|Fragment.ID), data = master)
# qqnorm(resid(TP_CS_lmer))
# qqline(resid(TP_CS_lmer))
# 
# boxplot(resid(TP_CS_lmer)~master$Group)
# boxplot(resid(TP_CS_lmer)~master$Day)
anova(TP_CS_lmer, type="II")
## Type II Analysis of Variance Table with Satterthwaite's method
##             Sum Sq  Mean Sq NumDF DenDF F value   Pr(>F)   
## Group     0.046434 0.023217     2    12  1.4402 0.275047   
## Day       0.281598 0.140799     2    24  8.7339 0.001413 **
## Group:Day 0.085401 0.021350     4    24  1.3244 0.289469   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#capture.output(anova(TP_CS_lmer, type="II"), file = "../output/Statistics/TP_CS_lmer.csv")

###############################

Sym_Box <- ggplot(master, aes(x=Day, y=Cells.cm2.x6, fill = Group)) +
  geom_boxplot(width=.5, outlier.shape= NA, position = position_dodge(width = 0.5)) +
  stat_summary(fun=mean, geom="line", aes(group=Group, color = Group), position = position_dodge(width = 0.5))  + 
  #  stat_summary(fun=mean, geom="point")
  geom_point(pch = 21, position=position_jitterdodge(dodge.width=0.5)) +
  #  ylim(0,0.5) +
  scale_fill_manual(values=c("#E7B800", "#00AFBB", "#FC4E07")) +
  scale_color_manual(values=c("#E7B800", "#00AFBB", "#FC4E07")) + 
  scale_x_discrete(labels=c("0" = "Day 0", "37" = "Day 37", "52" = "Day 52")) +
  xlab("Timepoint") +  ylab(expression("Cell Density " (10^{6} ~ cm^{-2})))+ #label y axis + #Axis titles
  theme_bw() + theme(panel.border = element_rect(color="black", fill=NA, size=0.75), panel.grid.major = element_blank(), #Makes background theme white
                     panel.grid.minor = element_blank(), axis.line = element_blank()) +
  theme(axis.text = element_text(size = 15, color = "black"),
        axis.title = element_text(size = 18, color = "black"))+ ggtitle("Endosymbiont Cell Density (x10^6)")
Sym_Box

Sym_lmer <- lmer(Cells.cm2.x6~Group*Day+(1|Fragment.ID), data = master)
# qqnorm(resid(Sym_lmer))
# qqline(resid(Sym_lmer))
# 
# boxplot(resid(Sym_lmer)~master$Group)
# boxplot(resid(Sym_lmer)~master$Day)
anova(Sym_lmer, type="II")
## Type II Analysis of Variance Table with Satterthwaite's method
##            Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
## Group      2.4141  1.2071     2    12  17.927 0.0002486 ***
## Day       10.9398  5.4699     2    24  81.240 2.065e-11 ***
## Group:Day  0.4764  0.1191     4    24   1.769 0.1680518    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#capture.output(anova(Sym_lmer, type="II"), file = "../output/Statistics/Sym_lmer.csv")

#################################

Chla_SA_Box <- ggplot(master, aes(x=Day, y=ChlA.ugcm2, fill = Group)) +
  geom_boxplot(width=.5, outlier.shape= NA, position = position_dodge(width = 0.5)) +
  stat_summary(fun=mean, geom="line", aes(group=Group, color = Group), position = position_dodge(width = 0.5))  + 
  #  stat_summary(fun=mean, geom="point")
  geom_point(pch = 21, position=position_jitterdodge(dodge.width=0.5)) +
  #  ylim(0,0.5) +
  scale_fill_manual(values=c("#E7B800", "#00AFBB", "#FC4E07")) +
  scale_color_manual(values=c("#E7B800", "#00AFBB", "#FC4E07")) + 
  scale_x_discrete(labels=c("0" = "Day 0", "37" = "Day 37", "52" = "Day 52")) +
  xlab("Timepoint") +  ylab(expression("Chlorophyll a " (ug ~ cm^{-2})))+ #label y axis + #Axis titles
  theme_bw() + theme(panel.border = element_rect(color="black", fill=NA, size=0.75), panel.grid.major = element_blank(), #Makes background theme white
                     panel.grid.minor = element_blank(), axis.line = element_blank()) +
  theme(axis.text = element_text(size = 15, color = "black"),
        axis.title = element_text(size = 18, color = "black")) + ggtitle("Chlorophyll a (per cm2)")
Chla_SA_Box

Chla_SA_lmer <- lmer(ChlA.ugcm2~Group*Day+(1|Fragment.ID), data = master)
# qqnorm(resid(Chla_SA_lmer))
# qqline(resid(Chla_SA_lmer)) # double check normality
# 
# boxplot(resid(Chla_SA_lmer)~master$Group)
# boxplot(resid(Chla_SA_lmer)~master$Day)
anova(Chla_SA_lmer, type="II")
## Type II Analysis of Variance Table with Satterthwaite's method
##            Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
## Group      49.211  24.606     2    36  6.5346  0.003792 ** 
## Day       288.137 144.068     2    36 38.2610 1.234e-09 ***
## Group:Day   2.627   0.657     4    36  0.1744  0.950087    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#capture.output(anova(Chla_SA_lmer, type="II"), file = "../output/Statistics/Chla_SA_lmer.csv")

#################################

ChlC2_SA_Box <- ggplot(master, aes(x=Day, y=ChlC2.ugcm2, fill = Group)) +
  geom_boxplot(width=.5, outlier.shape= NA, position = position_dodge(width = 0.5)) +
  stat_summary(fun=mean, geom="line", aes(group=Group, color = Group), position = position_dodge(width = 0.5))  + 
  #  stat_summary(fun=mean, geom="point")
  geom_point(pch = 21, position=position_jitterdodge(dodge.width=0.5)) +
  #  ylim(0,0.5) +
  scale_fill_manual(values=c("#E7B800", "#00AFBB", "#FC4E07")) +
  scale_color_manual(values=c("#E7B800", "#00AFBB", "#FC4E07")) + 
  scale_x_discrete(labels=c("0" = "Day 0", "37" = "Day 37", "52" = "Day 52")) +
  xlab("Timepoint") +  ylab(expression("Chlorophyll c2 " (ug ~ cm^{-2})))+ #label y axis + #Axis titles
  theme_bw() + theme(panel.border = element_rect(color="black", fill=NA, size=0.75), panel.grid.major = element_blank(), #Makes background theme white
                     panel.grid.minor = element_blank(), axis.line = element_blank()) +
  theme(axis.text = element_text(size = 15, color = "black"),
        axis.title = element_text(size = 18, color = "black")) + ggtitle("Chlorophyll c2 (per cm2)")
ChlC2_SA_Box

ChlC2_SA_lmer <- lmer(ChlC2.ugcm2~Group*Day+(1|Fragment.ID), data = master)
# qqnorm(resid(ChlC2_SA_lmer))
# qqline(resid(ChlC2_SA_lmer)) # double check normality
# 
# boxplot(resid(ChlC2_SA_lmer)~master$Group)
# boxplot(resid(ChlC2_SA_lmer)~master$Day)
anova(ChlC2_SA_lmer, type="II")
## Type II Analysis of Variance Table with Satterthwaite's method
##            Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
## Group      4.3301  2.1650     2    36  1.4932 0.2382414    
## Day       31.0503 15.5251     2    36 10.7074 0.0002244 ***
## Group:Day  5.7109  1.4277     4    36  0.9847 0.4281952    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#capture.output(anova(ChlC2_SA_lmer, type="II"), file = "../output/Statistics/ChlC2_SA_lmer.csv")

#################################

Chla_cell_Box <- ggplot(master, aes(x=Day, y=ChlA.ugcell, fill = Group)) +
  geom_boxplot(width=.5, outlier.shape= NA, position = position_dodge(width = 0.5)) +
  stat_summary(fun=mean, geom="line", aes(group=Group, color = Group), position = position_dodge(width = 0.5))  + 
  #  stat_summary(fun=mean, geom="point")
  geom_point(pch = 21, position=position_jitterdodge(dodge.width=0.5)) +
  #  ylim(0,0.5) +
  scale_fill_manual(values=c("#E7B800", "#00AFBB", "#FC4E07")) +
  scale_color_manual(values=c("#E7B800", "#00AFBB", "#FC4E07")) + 
  scale_x_discrete(labels=c("0" = "Day 0", "37" = "Day 37", "52" = "Day 52")) +
  xlab("Timepoint") +  ylab(expression("Chlorophyll a " (ug ~ cell^{-1})))+ #label y axis + #Axis titles
  theme_bw() + theme(panel.border = element_rect(color="black", fill=NA, size=0.75), panel.grid.major = element_blank(), #Makes background theme white
                     panel.grid.minor = element_blank(), axis.line = element_blank()) +
  theme(axis.text = element_text(size = 15, color = "black"),
        axis.title = element_text(size = 18, color = "black")) + ggtitle("Chlorophyll a (per cell)")
Chla_cell_Box

Chla_cell_lmer <- lmer(ChlA.ugcell~Group*Day+(1|Fragment.ID), data = master)
# qqnorm(resid(Chla_cell_lmer))
# qqline(resid(Chla_cell_lmer)) # double check normality
# 
# boxplot(resid(Chla_cell_lmer)~master$Group)
# boxplot(resid(Chla_cell_lmer)~master$Day)
anova(Chla_cell_lmer, type="II")
## Type II Analysis of Variance Table with Satterthwaite's method
##               Sum Sq    Mean Sq NumDF DenDF F value Pr(>F)
## Group     1.3966e-11 6.9832e-12     2 1.141  2.8305 0.3611
## Day       5.5439e-12 2.7719e-12     2 1.141  1.1235 0.5375
## Group:Day 1.7957e-11 4.4892e-12     4 1.141  1.8196 0.4774
#capture.output(anova(Chla_cell_lmer, type="II"), file = "../output/Statistics/Chla_cell_lmer.csv")

#################################

ChlC2_cell_Box <- ggplot(master, aes(x=Day, y=ChlC2.ugcell, fill = Group)) +
  geom_boxplot(width=.5, outlier.shape= NA, position = position_dodge(width = 0.5)) +
  stat_summary(fun=mean, geom="line", aes(group=Group, color = Group), position = position_dodge(width = 0.5))  + 
  #  stat_summary(fun=mean, geom="point")
  geom_point(pch = 21, position=position_jitterdodge(dodge.width=0.5)) +
  #  ylim(0,0.5) +
  scale_fill_manual(values=c("#E7B800", "#00AFBB", "#FC4E07")) +
  scale_color_manual(values=c("#E7B800", "#00AFBB", "#FC4E07")) + 
  scale_x_discrete(labels=c("0" = "Day 0", "37" = "Day 37", "52" = "Day 52")) +
  xlab("Timepoint") +  ylab(expression("Chlorophyll c2 " (ug ~ cell^{-1})))+ #label y axis + #Axis titles
  theme_bw() + theme(panel.border = element_rect(color="black", fill=NA, size=0.75), panel.grid.major = element_blank(), #Makes background theme white
                     panel.grid.minor = element_blank(), axis.line = element_blank()) +
  theme(axis.text = element_text(size = 15, color = "black"),
        axis.title = element_text(size = 18, color = "black")) + ggtitle("Chlorophyll C2 (per cell)")
ChlC2_cell_Box

ChlC2_cell_lmer <- lmer(ChlC2.ugcell~Group*Day+(1|Fragment.ID), data = master)
# qqnorm(resid(ChlC2_cell_lmer))
# qqline(resid(ChlC2_cell_lmer)) # double check normality
# 
# boxplot(resid(ChlC2_cell_lmer)~master$Group)
# boxplot(resid(ChlC2_cell_lmer)~master$Day)
anova(ChlC2_cell_lmer, type="II")
## Type II Analysis of Variance Table with Satterthwaite's method
##               Sum Sq    Mean Sq NumDF  DenDF F value Pr(>F)
## Group     6.9132e-12 3.4566e-12     2 1.0236  1.6444 0.4790
## Day       3.2920e-13 1.6460e-13     2 1.0236  0.0783 0.9297
## Group:Day 1.4039e-11 3.5096e-12     4 1.0236  1.6697 0.5139
#capture.output(anova(ChlC2_cell_lmer, type="II"), file = "../output/Statistics/ChlC2_cell_lmer.csv")

#################################

Carb_Coral_Box <- ggplot(master, aes(x=Day, y=Coral_Carb_mgcm2, fill = Group)) +
  geom_boxplot(width=.5, outlier.shape= NA, position = position_dodge(width = 0.5)) +
  stat_summary(fun=mean, geom="line", aes(group=Group, color = Group), position = position_dodge(width = 0.5))  + 
  #  stat_summary(fun=mean, geom="point")
  geom_point(pch = 21, position=position_jitterdodge(dodge.width=0.5)) +
  #  ylim(0,0.5) +
  scale_fill_manual(values=c("#E7B800", "#00AFBB", "#FC4E07")) +
  scale_color_manual(values=c("#E7B800", "#00AFBB", "#FC4E07")) + 
  scale_x_discrete(labels=c("0" = "Day 0", "37" = "Day 37", "52" = "Day 52")) +
  xlab("Timepoint") +  ylab(expression("Total Carbohydrate " (mg~cm^{-2}))) + #Axis titles
  theme_bw() + theme(panel.border = element_rect(color="black", fill=NA, size=0.75), panel.grid.major = element_blank(), #Makes background theme white
                     panel.grid.minor = element_blank(), axis.line = element_blank()) +
  theme(axis.text = element_text(size = 15, color = "black"),
        axis.title = element_text(size = 18, color = "black")) + ggtitle("Coral Total Carbohydrate (per cm2)")
Carb_Coral_Box

Carb_Coral_lmer <- lmer(Coral_Carb_mgcm2~Group*Day+(1|Fragment.ID), data = master)
# qqnorm(resid(Carb_Coral_lmer))
# qqline(resid(Carb_Coral_lmer)) # double check normality
# 
# boxplot(resid(Carb_Coral_lmer)~master$Group)
# boxplot(resid(Carb_Coral_lmer)~master$Day)
anova(Carb_Coral_lmer, type="II")
## Type II Analysis of Variance Table with Satterthwaite's method
##           Sum Sq Mean Sq NumDF  DenDF F value   Pr(>F)   
## Group     19.129   9.565     2 11.992  2.1777 0.156027   
## Day       71.143  35.572     2 23.992  8.0990 0.002053 **
## Group:Day 20.127   5.032     4 23.992  1.1456 0.359181   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#capture.output(anova(Carb_Coral_lmer, type="II"), file = "../output/Statistics/Carb_Coral_lmer.csv")

#################################

Carb_Sym_SA_Box <- ggplot(master, aes(x=Day, y=Symbiont_Carb_mgcm2, fill = Group)) +
  geom_boxplot(width=.5, outlier.shape= NA, position = position_dodge(width = 0.5)) +
  stat_summary(fun=mean, geom="line", aes(group=Group, color = Group), position = position_dodge(width = 0.5))  + 
  #  stat_summary(fun=mean, geom="point")
  geom_point(pch = 21, position=position_jitterdodge(dodge.width=0.5)) +
  #  ylim(0,0.5) +
  scale_fill_manual(values=c("#E7B800", "#00AFBB", "#FC4E07")) +
  scale_color_manual(values=c("#E7B800", "#00AFBB", "#FC4E07")) + 
  scale_x_discrete(labels=c("0" = "Day 0", "37" = "Day 37", "52" = "Day 52")) +
  xlab("Timepoint") +  ylab(expression("Total Carbohydrate " (mg~cm^{-2}))) + #Axis titles
  theme_bw() + theme(panel.border = element_rect(color="black", fill=NA, size=0.75), panel.grid.major = element_blank(), #Makes background theme white
                     panel.grid.minor = element_blank(), axis.line = element_blank()) +
  theme(axis.text = element_text(size = 15, color = "black"),
        axis.title = element_text(size = 18, color = "black"))  + ggtitle("Symbiont Total Carbohydrate (per cm2)")
Carb_Sym_SA_Box

Carb_Sym_SA_lmer <- lmer(Symbiont_Carb_mgcm2~Group*Day+(1|Fragment.ID), data = master)
# qqnorm(resid(Carb_Sym_SA_lmer))
# qqline(resid(Carb_Sym_SA_lmer)) # double check normality
# 
# boxplot(resid(Carb_Sym_SA_lmer)~master$Group)
# boxplot(resid(Carb_Sym_SA_lmer)~master$Day)
anova(Carb_Sym_SA_lmer, type="II")
## Type II Analysis of Variance Table with Satterthwaite's method
##            Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Group     0.30567 0.15283     2    12  0.1354 0.8747
## Day       1.87116 0.93558     2    24  0.8290 0.4486
## Group:Day 1.39687 0.34922     4    24  0.3094 0.8688
#capture.output(anova(Carb_Sym_SA_lmer, type="II"), file = "../output/Statistics/Carb_Sym_SA_lmer.csv")

#################################

Carb_Sym_Cell_Box <- ggplot(master, aes(x=Day, y=Carb.mgcell, fill = Group)) +
  geom_boxplot(width=.5, outlier.shape= NA, position = position_dodge(width = 0.5)) +
  stat_summary(fun=mean, geom="line", aes(group=Group, color = Group), position = position_dodge(width = 0.5))  + 
  #  stat_summary(fun=mean, geom="point")
  geom_point(pch = 21, position=position_jitterdodge(dodge.width=0.5)) +
  #  ylim(0,0.5) +
  scale_fill_manual(values=c("#E7B800", "#00AFBB", "#FC4E07")) +
  scale_color_manual(values=c("#E7B800", "#00AFBB", "#FC4E07")) + 
  scale_x_discrete(labels=c("0" = "Day 0", "37" = "Day 37", "52" = "Day 52")) +
  xlab("Timepoint") +  ylab(expression("Total Carbohydrate " (mg~cm^{-2}))) + #Axis titles
  theme_bw() + theme(panel.border = element_rect(color="black", fill=NA, size=0.75), panel.grid.major = element_blank(), #Makes background theme white
                     panel.grid.minor = element_blank(), axis.line = element_blank()) +
  theme(axis.text = element_text(size = 15, color = "black"),
        axis.title = element_text(size = 18, color = "black"))  + ggtitle("Symbiont Total Carbohydrate (per cell)")
Carb_Sym_Cell_Box

Carb_Sym_Cell_lmer <- lmer(Carb.mgcell~Group*Day+(1|Fragment.ID), data = master)
# qqnorm(resid(Carb_Sym_Cell_lmer))
# qqline(resid(Carb_Sym_Cell_lmer)) # double check normality
# 
# boxplot(resid(Carb_Sym_Cell_lmer)~master$Group)
# boxplot(resid(Carb_Sym_Cell_lmer)~master$Day)
anova(Carb_Sym_Cell_lmer, type="II")
## Type II Analysis of Variance Table with Satterthwaite's method
##               Sum Sq    Mean Sq NumDF  DenDF F value  Pr(>F)  
## Group     3.2244e-10 1.6122e-10     2 2.2679 12.7207 0.05853 .
## Day       4.4263e-10 2.2131e-10     2 2.7318 17.4623 0.02778 *
## Group:Day 1.5965e-10 3.9911e-11     4 2.7318  3.1491 0.20124  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#capture.output(anova(Carb_Sym_Cell_lmer, type="II"), file = "../output/Statistics/Carb_Sym_Cell_lmer.csv")

#################################

Carb_CS_Box <- ggplot(master, aes(x=Day, y=HostSymRatioCarb, fill = Group)) +
  geom_boxplot(width=.5, outlier.shape= NA, position = position_dodge(width = 0.5)) +
  stat_summary(fun=mean, geom="line", aes(group=Group, color = Group), position = position_dodge(width = 0.5))  + 
  #  stat_summary(fun=mean, geom="point")
  geom_point(pch = 21, position=position_jitterdodge(dodge.width=0.5)) +
  #  ylim(0,0.5) +
  scale_fill_manual(values=c("#E7B800", "#00AFBB", "#FC4E07")) +
  scale_color_manual(values=c("#E7B800", "#00AFBB", "#FC4E07")) + 
  scale_x_discrete(labels=c("0" = "Day 0", "37" = "Day 37", "52" = "Day 52")) +
  xlab("Timepoint") +  ylab("Coral to Symbiont Ratio (Total Carbohydrate)")  + #Axis titles
  theme_bw() + theme(panel.border = element_rect(color="black", fill=NA, size=0.75), panel.grid.major = element_blank(), #Makes background theme white
                     panel.grid.minor = element_blank(), axis.line = element_blank()) +
  theme(axis.text = element_text(size = 15, color = "black"),
        axis.title = element_text(size = 18, color = "black")) + ggtitle("Host to Symbiont Total Carbohydrate")
Carb_CS_Box

Carb_CS_lmer <- lmer(HostSymRatioCarb~Group*Day+(1|Fragment.ID), data = master)
# qqnorm(resid(Carb_CS_lmer))
# qqline(resid(Carb_CS_lmer)) # double check normality
# 
# boxplot(resid(Carb_CS_lmer)~master$Group)
# boxplot(resid(Carb_CS_lmer)~master$Day)
anova(Carb_CS_lmer, type="II")
## Type II Analysis of Variance Table with Satterthwaite's method
##           Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Group     12.947  6.4738     2    12  1.3349 0.2996
## Day       10.740  5.3701     2    24  1.1073 0.3468
## Group:Day 10.730  2.6824     4    24  0.5531 0.6987
#capture.output(anova(Carb_CS_lmer, type="II"), file = "../output/Statistics/Carb_CS_lmer.csv")

#################################

PG_Box <- ggplot(master, aes(x=Day, y=Pgross_umol.cm2.hr, fill = Group)) +
  geom_boxplot(width=.5, outlier.shape= NA, position = position_dodge(width = 0.5)) +
  stat_summary(fun=mean, geom="line", aes(group=Group, color = Group), position = position_dodge(width = 0.5))  + 
#  stat_summary(fun=mean, geom="point")
  geom_point(pch = 21, position=position_jitterdodge(dodge.width=0.5), outlier.shape= NA) +
#  ylim(0,0.5) +
  scale_fill_manual(values=c("#E7B800", "#00AFBB", "#FC4E07")) +
  scale_color_manual(values=c("#E7B800", "#00AFBB", "#FC4E07")) + 
  scale_x_discrete(labels=c("T1" = "Day 0", "T2" = "Day 37", "T3" = "Day 52")) +
  xlab("Timepoint") +  ylab(expression("Gross Photosynthetic Rate " (mu*mol ~ cm^{-2} ~ h^{-1}))) + #label y axis + #Axis titles
  theme_bw() + theme(panel.border = element_rect(color="black", fill=NA, size=0.75), panel.grid.major = element_blank(), #Makes background theme white
                     panel.grid.minor = element_blank(), axis.line = element_blank()) +
  theme(axis.text = element_text(size = 15, color = "black"),
        axis.title = element_text(size = 18, color = "black")) + ggtitle("Gross Photosynthetic Rate")
PG_Box

PG_lmer <- lmer(Pgross_umol.cm2.hr~Group*Day+(1|Fragment.ID), data = master)
# qqnorm(resid(PG_lmer))
# qqline(resid(PG_lmer)) # double check normality
# 
# boxplot(resid(PG_lmer)~master$Group)
# boxplot(resid(PG_lmer)~master$Day)
anova(PG_lmer, type="II")
## Type II Analysis of Variance Table with Satterthwaite's method
##           Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
## Group     0.6182 0.30910     2    12  14.942  0.000553 ***
## Day       3.2558 1.62790     2    24  78.696 2.878e-11 ***
## Group:Day 1.3409 0.33522     4    24  16.205 1.478e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#capture.output(anova(PG_lmer, type="II"), file = "../output/Statistics/PG_lmer.csv")

#################################

Resp_Box <- ggplot(master, aes(x=Day, y=abs.Rdark_umol.cm2.hr, fill = Group)) +
  geom_boxplot(width=.5, outlier.shape= NA, position = position_dodge(width = 0.5)) +
  stat_summary(fun=mean, geom="line", aes(group=Group, color = Group), position = position_dodge(width = 0.5))  + 
  #  stat_summary(fun=mean, geom="point")
  geom_point(pch = 21, position=position_jitterdodge(dodge.width=0.5), outlier.shape= NA) +
  #  ylim(0,0.5) +
  scale_fill_manual(values=c("#E7B800", "#00AFBB", "#FC4E07")) +
  scale_color_manual(values=c("#E7B800", "#00AFBB", "#FC4E07")) + 
  scale_x_discrete(labels=c("T1" = "Day 0", "T2" = "Day 37", "T3" = "Day 52")) +
  xlab("Timepoint") +  ylab(expression("Respiration Rate " (mu*mol ~ cm^{-2} ~ h^{-1}))) + #label y axis + #Axis titles
  theme_bw() + theme(panel.border = element_rect(color="black", fill=NA, size=0.75), panel.grid.major = element_blank(), #Makes background theme white
                     panel.grid.minor = element_blank(), axis.line = element_blank()) +
  theme(axis.text = element_text(size = 15, color = "black"),
        axis.title = element_text(size = 18, color = "black")) + ggtitle("Light-Enhanced Respiration Rate")
Resp_Box

Resp_lmer <- lmer(abs.Rdark_umol.cm2.hr~Group*Day+(1|Fragment.ID), data = master)
# qqnorm(resid(Resp_lmer))
# qqline(resid(Resp_lmer)) # double check normality
# 
# boxplot(resid(Resp_lmer)~master$Group)
# boxplot(resid(Resp_lmer)~master$Day)
anova(Resp_lmer, type="II")
## Type II Analysis of Variance Table with Satterthwaite's method
##            Sum Sq  Mean Sq NumDF DenDF F value    Pr(>F)    
## Group     0.04322 0.021610     2    36  1.9296    0.1599    
## Day       0.57080 0.285398     2    36 25.4839 1.273e-07 ***
## Group:Day 0.03669 0.009172     4    36  0.8190    0.5216    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#capture.output(anova(Resp_lmer, type="II"), file = "../output/Statistics/Resp_lmer.csv")

#################################

PR_Box <- ggplot(master, aes(x=Day, y=PR, fill = Group)) +
  geom_boxplot(width=.5, outlier.shape= NA, position = position_dodge(width = 0.5)) +
  stat_summary(fun=mean, geom="line", aes(group=Group, color = Group), position = position_dodge(width = 0.5))  + 
  #  stat_summary(fun=mean, geom="point")
  geom_point(pch = 21, position=position_jitterdodge(dodge.width=0.5), outlier.shape= NA) +
  #  ylim(0,0.5) +
  scale_fill_manual(values=c("#E7B800", "#00AFBB", "#FC4E07")) +
  scale_color_manual(values=c("#E7B800", "#00AFBB", "#FC4E07")) + 
  scale_x_discrete(labels=c("T1" = "Day 0", "T2" = "Day 37", "T3" = "Day 52")) +
  xlab("Timepoint") +  ylab("P:R") + #label y axis + #Axis titles
  theme_bw() + theme(panel.border = element_rect(color="black", fill=NA, size=0.75), panel.grid.major = element_blank(), #Makes background theme white
                     panel.grid.minor = element_blank(), axis.line = element_blank()) +
  theme(axis.text = element_text(size = 15, color = "black"),
        axis.title = element_text(size = 18, color = "black")) + ggtitle("Photosyntheis to Respiration Ratio")
PR_Box

PR_lmer <- lmer(PR~Group*Day+(1|Fragment.ID), data = master)
# qqnorm(resid(PR_lmer))
# qqline(resid(PR_lmer)) # double check normality
# 
# boxplot(resid(PR_lmer)~master$Group)
# boxplot(resid(PR_lmer)~master$Day)
anova(PR_lmer, type="II")
## Type II Analysis of Variance Table with Satterthwaite's method
##           Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
## Group     4.3386  2.1693     2    12  41.894 3.866e-06 ***
## Day       4.2652  2.1326     2    24  41.185 1.741e-08 ***
## Group:Day 4.3340  1.0835     4    24  20.924 1.549e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#capture.output(anova(PR_lmer, type="II"), file = "../output/Statistics/PR_lmer.csv")

#################################

PCA

master.pca <- master %>%
  dplyr::select(-(Surface_Area_cm2)) 

scaled_pca <-prcomp(master.pca[c(4:20)], scale=TRUE, center=TRUE)
fviz_eig(scaled_pca)

coral_info<-master.pca[c(2,3)]

pca_data<- scaled_pca%>%
  augment(coral_info)%>%
  group_by(Day, Group)%>%
  mutate(PC1.mean = mean(.fittedPC1),
         PC2.mean = mean(.fittedPC2))

pca.centroids<- pca_data%>% 
  dplyr::select(Day, Group, PC1.mean, PC2.mean)%>%
  dplyr::group_by(Day, Group)%>%
  dplyr::summarise(PC1.mean = mean(PC1.mean),
                   PC2.mean = mean(PC2.mean))

#Examine PERMANOVA results.  

# scale data
vegan <- scale(master.pca[c(4:20)])

# PerMANOVA 
permanova<-adonis(vegan ~ Day*Group, data = master.pca, method='eu')
z_pca<-permanova$aov.tab
z_pca
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## Day        2    200.80 100.400  9.4360 0.26845  0.001 ***
## Group      2     98.64  49.322  4.6354 0.13188  0.001 ***
## Day:Group  4     65.51  16.377  1.5392 0.08758  0.074 .  
## Residuals 36    383.05  10.640         0.51209           
## Total     44    748.00                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Assemble plot with background points

#1. make plot with dots

#adding percentages on axis

names(pca_data)[4] <- "PCA1"
names(pca_data)[5] <- "PCA2"
percentage <- round((scaled_pca$sdev^2) / sum((scaled_pca$sdev^2)) * 100, 2)
percentage <- paste( colnames(pca_data[4:20]), "(", paste(as.character(percentage), "%", ")", sep="") )

#setting up data to add polygons
pca_data$Day.Group <- paste(pca_data$Day, pca_data$Group)
find_hull <- function(pca_data) pca_data[chull(pca_data$PCA1, pca_data$PCA2), ]
hulls <- ddply(pca_data, "Day.Group", find_hull)

PCA<-ggplot(pca_data, aes(PCA1, PCA2, color=Group)) + 
  geom_point(size = 4, alpha=0.2, aes(shape = Day))+
  scale_colour_manual(values=c("#46008B", "#8B0046", "#468B00")) +
  scale_fill_manual(values=c("#46008B", "#8B0046", "#468B00")) + 
  scale_shape_manual(values=c(15, 17, 19)) +
  theme_classic()+
  ylim(-7,12)+
  xlim(-6,6)+
  ylab(percentage[2])+
  xlab(percentage[1])+
  geom_text(x=4.5, y=-4.75, label=paste("p(Day)=", z_pca$`Pr(>F)`[1]), size=4, color=ifelse(z_pca$`Pr(>F)`[1] < 0.05, "black", "darkgray")) + 
  geom_text(x=4.5, y=-5.5, label=paste("p(Group)=", z_pca$`Pr(>F)`[2]), size=4, color=ifelse(z_pca$`Pr(>F)`[2] < 0.05, "black", "darkgray")) + 
  geom_text(x=4.5, y=-6.25, label=paste("p(Day x Group)=", z_pca$`Pr(>F)`[3]), size=4, color=ifelse(z_pca$`Pr(>F)`[3] < 0.05, "black", "darkgray")) + 
  theme(legend.text = element_text(size=8), 
        legend.position="none",
        plot.background = element_blank(),
        legend.title = element_text(size=10), 
        plot.margin = margin(1, 1, 1, 1, "cm"),
        axis.text = element_text(size=18), 
        title = element_text(size=25, face="bold"), 
        axis.title = element_text(size=18))

#Add centroids  

#2. add centroids 
PCAcen<-PCA +  geom_polygon(data = hulls, alpha = 0.2, aes(color = Group, fill = Group, lty = Day)) +
  geom_point(aes(x=PC1.mean, y=PC2.mean,color=Group, shape = Day), data=pca.centroids, size=4, show.legend=FALSE) + 
  scale_linetype_manual(values = c("solid", "dashed", "dotted")) +
  scale_colour_manual(values=c("#46008B", "#8B0046", "#468B00"), breaks = c("Control","Bleached", "Mortality"), labels = c("Control", "Bleached", "Partial Mortality")) +
  scale_fill_manual(values=c("#46008B", "#8B0046", "#468B00"), breaks = c("Control","Bleached", "Mortality"), labels = c("Control", "Bleached", "Partial Mortality")) + 
  scale_shape_manual(values=c(15, 17, 19)) +
  theme(legend.text = element_text(size=8), 
        legend.position=c(0.95,0.85),
        plot.background = element_blank(),
        legend.title = element_text(size=10), 
        plot.margin = margin(1, 1, 1, 1, "cm"),
        axis.text = element_text(size=18), 
        title = element_text(size=25, face="bold"), 
        axis.title = element_text(size=18))

#Add segments

#3. add segments
segpoints<-pca.centroids%>%
  gather(variable, value, -(Day:Group)) %>%
  unite(temp, Day, variable) %>%
  spread(temp, value) 
  
names(segpoints)[2] <- "Day0_PC1.mean"
names(segpoints)[3] <- "Day0_PC2.mean"
names(segpoints)[4] <- "Day37_PC1.mean"
names(segpoints)[5] <- "Day37_PC2.mean"
names(segpoints)[6] <- "Day52_PC1.mean"
names(segpoints)[7] <- "Day52_PC2.mean"

PCAfull<-PCAcen + 
  geom_segment(aes(x = Day0_PC1.mean, y = Day0_PC2.mean, xend = Day37_PC1.mean, yend = Day37_PC2.mean, colour = Group), data = segpoints, size=2, show.legend=FALSE) +
  geom_segment(aes(x = Day37_PC1.mean, y = Day37_PC2.mean, xend = Day52_PC1.mean, yend = Day52_PC2.mean, colour = Group), data = segpoints, size=2, arrow = arrow(length=unit(0.5,"cm")), show.legend=FALSE); PCAfull

#Add bi plot with loadings  

#1. make plot with dots
biplotArrows<-ggplot2::autoplot(scaled_pca, data=master.pca, loadings=TRUE,  colour="Group", loadings.label.colour="black", loadings.colour="black", loadings.label=TRUE, frame=FALSE, loadings.label.size=3.5, loadings.label.vjust=0.5, loadings.label.hjust=-0.1, loadings.label.repel=FALSE, size=4, alpha=0.0) + 
  #scale_colour_manual(values=c("blue4", "darkgray", "springgreen4")) +
  #scale_fill_manual(values=c("blue4", "darkgray", "springgreen4")) + 
  theme_classic()+
  #xlim(-.3,.6)+
  #ylim(-.3, .3)+
  theme(legend.text = element_text(size=18), 
        legend.position="none",
        axis.text.x=element_blank(),
        axis.text.y=element_blank(),
        plot.background = element_blank(),
        legend.title = element_blank(), 
        plot.margin = margin(1, 1, 1, 1, "cm"),
        axis.text = element_text(size=18), 
        title = element_text(size=25, face="bold"),
        axis.title = element_blank());biplotArrows

FinalFullPCA<-ggdraw(PCAfull) + #theme_half_open(12)) +
  draw_plot(biplotArrows, .05, .6, .4, .4); FinalFullPCA #x, y, w, h

ggsave(filename="../output/Physiology/FullPCA_phys.pdf", plot=FinalFullPCA, dpi=300, width=12, height=10, units="in")